% COZZI AND IMPULLITTI (2015), GLOBALIZATION AND WAGE POLARIZATION, REVIEW
% OF ECONOMICS AND STATISTICS

% CALIBRATION FUNCTION FILE 1: EQUILIBRIUM SYSTEM


function f=funz(x)


global sig1 gam beta phi AD AF epsi fi1 n lambda eta rho eserv fi tau zf1 zd1

the0D = x(1);
the0F = x(2);
wDL = x(3);

wDH = (sig1*(the0D/(the0D-gam))*x(3));
wFL = (((the0F-gam)/(sig1*the0F))^(1-beta));
wFH = (sig1*(the0F/(the0F-gam))*wFL);


% =========================== TECHNOLOGY  =================================

ZD = ((1/zd1)*(wDL^beta)*(wDH^(1-beta))); % Production cost function D
ZF = ((1/zf1)*(wFL^beta)*(wFH^(1-beta))); % Production cost function F
FD = ((1/AD)*((wDL^phi)*(wDH)^(1-phi)));  % Innovation cost function D
FF = ((1/AF)*((wFL^phi)*(wFH)^(1-phi)));  % Innovation cost function F

   % Marginal products: production
   
ZDL = (1/zd1)*(beta*((wDH/wDL)^(1-beta)));
ZDH = (1/zd1)*((1-beta)*((wDH/wDL)^(-beta)));
ZFL = (1/zf1)*(beta*((wFH/wFL)^(1-beta)));
ZFH = (1/zf1)*((1-beta)*((wFH/wFL)^(-beta)));
   
   % Marginal products: R&D
   
FDL = (phi*(1/AD)*((wDH/wDL)^(1-phi)));
FDH = ((1-phi)*(1/AD)*((wDH/wDL)^(-phi)));
FFL = (phi*(1/AF)*((wFH/wFL)^(1-phi)));
FFH = ((1-phi)*(1/AF)*((wFH/wFL)^(-phi)));

% ============================== CUTOFFS ==================================

   % Skill cutoff to hire service workers
   
theHSD_fun = @(y) 1-y^epsi-(((wDH/wDL)*(1-eserv)*(y-gam))^epsi);
theHSF_fun = @(y) 1-y^epsi-(((wFH/wFL)*(1-eserv)*(y-gam))^epsi);

theHSD0 = 0.6;
theHSF0 = 0.6;

theHSD = csolve(theHSD_fun,theHSD0,[],1e-7,300);
theHSF = csolve(theHSF_fun,theHSF0,[],1e-7,300);

   % Skill cutoff for unskilled in production

theLSD = ((1-(theHSD^epsi))^(1/epsi));
theLSF = ((1-(theHSF^epsi))^(1/epsi));
  
   % Labor supplies
   
LDm = ((epsi/(1+epsi))*(the0D^(1+epsi)-theLSD^(1+epsi)));  % unskilled supply of production workers
LFm = ((epsi/(1+epsi))*(the0F^(1+epsi)-theLSF^(1+epsi)));


HD = ((1-the0D^epsi)*fi*(((epsi*((1-theHSD^(1+epsi))*(2 - eserv) + theHSD^(1+epsi) - the0D^(1+epsi)))/((1+epsi)*(1-the0D^epsi)))-gam)); % Skilled labor supply
HF = ((1-the0F^epsi)*fi*(((epsi*((1-theHSF^(1+epsi))*(2 - eserv) + theHSF^(1+epsi) - the0F^(1+epsi)))/((1+epsi)*(1-the0F^epsi)))-gam));

   % Average workers skills
   
theDm_avg = ((epsi/(1+epsi)) *((the0D)^(1+epsi)-(theLSD)^(1+epsi))/(the0D^epsi-theLSD^epsi)); % production
theFm_avg = ((epsi/(1+epsi)) *((the0F)^(1+epsi)-(theLSF)^(1+epsi))/(the0F^epsi-theLSF^epsi));


theDS_avg = ((epsi/(1+epsi)) *theLSD);  % Unskilled (services)
theFS_avg = ((epsi/(1+epsi)) *theLSF);


theDH_avg = (((epsi*((1-theHSD^(1+epsi)) + theHSD^(1+epsi) - the0D^(1+epsi)))/((1+epsi)*(1-the0D^epsi)))-gam);  % This average corrects for the increased supply of labour by skills: it is the average over the individuals of the averages per unit time, which means that those who work 2-eserv hours have their wages divided by (2-eserv). This is why we do not see (2-eserv) there, which of course we see it in the aggregate skill supply
theFH_avg = (((epsi*((1-theHSF^(1+epsi)) + theHSF^(1+epsi) - the0F^(1+epsi)))/((1+epsi)*(1-the0F^epsi)))-gam);  % This average corrects for the increased supply of labour by skills: it is the average over the individuals of the averages per unit time, which means that those who work 2-eserv hours have their wages divided by (2-eserv). This is why we do not see (2-eserv) there, which of course we see it in the aggregate skill supply

theDunskill_avg = ((epsi/(1+epsi))*the0D); % Average unskilled (service+production)

% ==================== OTHER ENDOGENOUS VARIABLES =========================

% Here we report some expressions we use the reduced the original
% equilibrium system to 3 equations.

qF_til = ((1-fi1)*((lambda^(eta-1))-1)/n);
Fr = (FD/FF);
cD_til =( (Fr*(lambda-1)-(lambda-ZD*tau)) / ((lambda-ZD/tau)-(lambda-1)*Fr) );
x_til = ((1/FF)* (((lambda-1)/lambda)*(1+cD_til))/ (rho+(n/((1-fi1)*((lambda^(eta-1))-1))) + (n*fi1/(1-fi1))));
IF = (1/qF_til)*(1/(cD_til+1));
CF = ((1/IF)*(LFm / (qF_til*(ZFL/lambda)*(1+cD_til) + x_til*FFL)));
qF = (qF_til*IF);
qD = (1-(qF));
ID = ((1/qF_til)-IF);
CD = (qD/qF)*CF;
xindex = ((x_til)*CF); % R&D Difficulty index

% ===================== EQUILIBRIUM SYSTEM ================================

f(1) = HD- qD*(ZDH/lambda)*(CD/tau + CF*tau)-ID*xindex*FDH; % market clearing skilled D
f(2) = LDm- qD*(ZDL/lambda)*(CD/tau +CF*tau)-ID*xindex*FDL; % market clearing unkilled D
f(3) = HF - (LFm*(ZFH*qF_til*(1+cD_til)/lambda + x_til*FFH) / ( qF_til*((ZFL/lambda)*(1+cD_til))+x_til*FFL ));



